function[y0,C]=Spline_interp(X,Y,s0，sN,x0)
%X，Y是已知插值点坐标
%s0，sN是两端的一次导数值
% x0是插值点
% y0是三次样条函数在×0处的值
%C是分段三次样条函数的系数
N=length(X);
C=zeros(4,N-1); h=zeros(1,N-1);
mu-zeros(1,N-1); Imt=zeros(1,N-1);
d=zcros(1,N); %d表示右端函数值
h=X(1.2:N)-X(1,1:N-1);
mu(1,N 1)=1;lmt(1,1)=1;
mu(1,1;N 2)=h(1,1:N-2)./(b(1,1:N-2)+h(1,2:N-1)):lmt(1,2:N-1)=h(1,2:N-1)./(h(1,1,N-2)+h(1,2:-1));
d(1，1)=6*((Y(1.2)-Y(1,1))/h(1,1)-s0)/h(1,1);
d(1,N)=6(sN (Y(1,N)-Y(1,N-1))/h(1,N-1))/h(1,N-1);
d(1，2:N-1)-6*((Y(1,3N)Y(1,2,N-1))./h(1，2：N-1)-...
(Y(1，2：N-1)-Y(1,1:N-2))./h(1,1:N 2))./(h(1,1:N-2)+...
+h(1,2:N-1));
bit=zeros(1.N-1);
bit(1,1)=lmt(1.1)/2;
for i=2:N-1
bit(1,i)=lmt(1,i)/(2-mu(1,i-1)bit(li-1));
end
y=zeros(1.N);
y(1,1)=d(1,1)/2;
for i=2:N
y(1.i)=(d(l,i)- mu(1,i-1)*y(1,i-1))/(2-mu(1,i-1)bit(1,i-1));
end
x=zeros(1.N);
x(1,N)=y(1,N);
for i=-1:-1:1
x(1,i)=y(1,i)-bit(l,i)*x(l,i+1);
end
v-zeros(1,N-1)；
v(1，1：N-1)=(Y(1，2：N)-Y(1，1；N-1))./h(1，1:N D；
C(4，：)-Y(1，1:N-1)；
C(3，：)=v-h，x(2xx(1，1；N-1)|x(1，2；N))/6；
C(2，：)=x(1，1：N-1)/2；
C(1，：)=(x(1，2：N)-x(1，1：N-1))./(6*h)；
if nargin＜5
y0-0；
else
for i-1:N-1
if xo>-X(1.j)8x0=X(1，i+1)
omg-xū-X(1，j)：
y0=((C(4，j)xomg+C(3，j))xomg-C(2，j))omg-C(1，j)；
cnd
end
end
